﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Nature;
using Nature.Numerics;
using System.Diagnostics;
using Nature.Numerics.Basic;

namespace Nature.Toolkit.GEqu
{
    public class HPIdealGasChemicalEquilibriumState : IChemicalEquilibriumState, IFunction
    {
        #region Private Fields

        double m_h0;
        double? m_p0;
        double? m_t, m_t0;
        IChemicalEquilibriumState m_gibbs;
        Newton<HPIdealGasChemicalEquilibriumState> m_newton;
        double[] m_x0;
        double[] m_hk;
        bool m_compositionInitialized = false;
        bool m_equilibrated = false;

        #endregion

        #region Private Methods

        void Equilibrate()
        {
            if (m_equilibrated) { return; }
            if (false == m_t0.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullTemperatureException(); }
            if (false == m_p0.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullPressureException(); }
            if (false == m_compositionInitialized) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullCompositionException(); }

            this.Model.IdealGasSpeciesThermoPotentials.MolarH(m_t0.Value, m_hk);
            double wsum = this.Model.CalcMixtureMolarMassFromX(m_x0);
            double hsum = 0.0;
            for (int isp = 0; isp < Nsp; ++isp) { hsum += m_x0[isp] * m_hk[isp]; }
            m_h0 = hsum / wsum;

            m_t = m_newton.Solve(m_t ?? m_t0.Value);
            m_equilibrated = true;
        }
        #endregion

        #region Constructors

        public HPIdealGasChemicalEquilibriumState(ModelObject model)
        {
            m_gibbs = new TPIdealGasChemicalEquilibriumState(model);
            this.Model = model;
            m_newton = new Newton<HPIdealGasChemicalEquilibriumState>(this);
            m_newton.AbsRadius = 100.0;
            m_newton.RelRadius = 0.1;
            m_x0 = new double[model.Nsp];
            m_hk = new double[model.Nsp];
        }

        #endregion

        #region Public Properties

        public ModelObject Model { get; private set; }

        public int Nsp { get { return this.Model.Nsp; } }

        #endregion

        #region IChEqCalculator Members

        double IChemicalEquilibriumState.T
        {
            get
            {
                Equilibrate();
                return m_t.Value;
            }
            set
            {
                if (m_t0 != value)
                {
                    m_equilibrated = false;
                    m_t0 = value;
                }
            }
        }

        double IChemicalEquilibriumState.P
        {
            get
            {
                if (m_p0.HasValue) { return m_p0.Value; }
                throw ChemicalEquilibriumStateNotInitializedException.CreateNullPressureException();
            }
            set
            {
                if (m_p0 != value)
                {
                    m_equilibrated = false;
                    m_p0 = value;
                    m_gibbs.P = value;
                }
            }
        }

        ReadOnlyArray<double> IChemicalEquilibriumState.X
        {
            get
            {
                Equilibrate();
                return m_gibbs.X;
            }
            set
            {
                m_compositionInitialized = true;
                
                bool changing = false;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    changing |= (m_x0[isp] != value[isp]);
                    m_x0[isp] = value[isp];                    
                }
                m_gibbs.X = value;
                if (changing)
                {
                    m_equilibrated = false;                    
                }
                 
            }
        }

        void IChemicalEquilibriumState.DropState()
        {
            m_t = null;
            m_equilibrated = false;
            m_t0 = null;
            m_p0 = null;
            m_compositionInitialized = false;
            m_gibbs.DropState();
        }

        #endregion

        #region IFunction Members

        double IFunction.Compute(double temperature)
        {
            m_gibbs.T = temperature;
            var x = m_gibbs.X;
            double wsum = this.Model.CalcMixtureMolarMassFromX(x);
            this.Model.IdealGasSpeciesThermoPotentials.MolarH(temperature, m_hk);
            double hsum = 0.0;
            for (int isp = 0; isp < Nsp; ++isp) { hsum += x[isp] * m_hk[isp]; }
            double diff = m_h0 - hsum / wsum;
            return diff;
        }

        #endregion
    }
}